library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(knitr)
library(ggplot2)
library(dplyr)
library(lubridate)
library(MMWRweek)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
library(readxl)
library(excessmort)
library(plotly)
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
library(patchwork)
library(scales)
## 
## Attaching package: 'scales'
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## The following object is masked from 'package:readr':
## 
##     col_factor
library(cowplot)
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:patchwork':
## 
##     align_plots
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
library(broom)
library(gt)
## Warning: Can't find generic `as.gtable` in package gtable to register S3 method.
## ℹ This message is only shown to developers using devtools.
## ℹ Do you need to update gtable to the latest version?
## 
## Attaching package: 'gt'
## 
## The following object is masked from 'package:cowplot':
## 
##     as_gtable
library(webshot2)
library(stringr)
library(purrr)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
source("~/Desktop/Thesis Work/Hospitalization/11-26 data/Deficit_Model_1.0.R")
allcause_state_byweek <- read_excel("allcause_state_byweek_suppressed_25mar21.xlsx")
icd10categories <- read_excel("icd10categories_updated_321.xlsx")

#weekly state data
all_cause_state_byweek <- allcause_state_byweek %>%
  rename(weekly_count_statelevel = count)

all_cause_state_by_week <- all_cause_state_byweek %>%
  separate(mmwr_week, into = c("MMWRyear", "MMWRweek"), sep = "-", convert = TRUE) %>%
  mutate(
    date = MMWRweek2Date(MMWRyear = MMWRyear, MMWRweek = MMWRweek, MMWRday = 1),
    statecount_numeric = as.numeric(ifelse(weekly_count_statelevel == "<5", NA, weekly_count_statelevel)),
    statecount_suppressed = weekly_count_statelevel == "<5",
    category = as.character(category)
  )

all_cause_state_by_week %>%
  group_by(category) %>%
  reframe(missing_total = sum(statecount_suppressed), 
            count_total = sum(statecount_numeric, na.rm = TRUE))
## # A tibble: 10 × 3
##    category               missing_total count_total
##    <chr>                          <int>       <dbl>
##  1 cardiovascular                     0      373147
##  2 covid                              1       58262
##  3 digestive                          0      182701
##  4 infectious                         0      112279
##  5 infectious respiratory             0       29816
##  6 injury                             0      331205
##  7 other                              0     1740863
##  8 pneumonia                          0       97954
##  9 renal                              0      102199
## 10 respiratory                        0      164424
all_cause_state_by_week
## # A tibble: 3,140 × 7
##    MMWRyear MMWRweek category               weekly_count_statelevel date      
##       <int>    <int> <chr>                  <chr>                   <date>    
##  1     2019        1 cardiovascular         843                     2018-12-30
##  2     2019        1 digestive              441                     2018-12-30
##  3     2019        1 infectious             299                     2018-12-30
##  4     2019        1 infectious respiratory 240                     2018-12-30
##  5     2019        1 injury                 691                     2018-12-30
##  6     2019        1 other                  4090                    2018-12-30
##  7     2019        1 pneumonia              427                     2018-12-30
##  8     2019        1 renal                  220                     2018-12-30
##  9     2019        1 respiratory            522                     2018-12-30
## 10     2019        1 covid                  0                       2018-12-30
## # ℹ 3,130 more rows
## # ℹ 2 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>
str(icd10categories)
## tibble [869 × 5] (S3: tbl_df/tbl/data.frame)
##  $ icd10code          : chr [1:869] "A00" "A01" "A02" "A03" ...
##  $ category_original  : chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
##  $ category_1113update: chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
##  $ category_123update : chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
##  $ category_0321update: chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
names(icd10categories)
## [1] "icd10code"           "category_original"   "category_1113update"
## [4] "category_123update"  "category_0321update"
dim(icd10categories)
## [1] 869   5
summary(icd10categories)
##   icd10code         category_original  category_1113update category_123update
##  Length:869         Length:869         Length:869          Length:869        
##  Class :character   Class :character   Class :character    Class :character  
##  Mode  :character   Mode  :character   Mode  :character    Mode  :character  
##  category_0321update
##  Length:869         
##  Class :character   
##  Mode  :character
sapply(icd10categories, class)
##           icd10code   category_original category_1113update  category_123update 
##         "character"         "character"         "character"         "character" 
## category_0321update 
##         "character"
glimpse(icd10categories)
## Rows: 869
## Columns: 5
## $ icd10code           <chr> "A00", "A01", "A02", "A03", "A04", "A05", "A06", "…
## $ category_original   <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_1113update <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_123update  <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_0321update <chr> "infectious", "infectious", "infectious", "infecti…
colSums(is.na(icd10categories))
##           icd10code   category_original category_1113update  category_123update 
##                   0                 418                 418                   1 
## category_0321update 
##                   0
str(allcause_state_byweek)
## tibble [3,140 × 3] (S3: tbl_df/tbl/data.frame)
##  $ mmwr_week: chr [1:3140] "2019-01" "2019-01" "2019-01" "2019-01" ...
##  $ category : chr [1:3140] "cardiovascular" "digestive" "infectious" "infectious respiratory" ...
##  $ count    : chr [1:3140] "843" "441" "299" "240" ...
names(allcause_state_byweek)
## [1] "mmwr_week" "category"  "count"
dim(allcause_state_byweek)
## [1] 3140    3
summary(allcause_state_byweek)
##   mmwr_week           category            count          
##  Length:3140        Length:3140        Length:3140       
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character
sapply(allcause_state_byweek,class)
##   mmwr_week    category       count 
## "character" "character" "character"
glimpse(allcause_state_byweek)
## Rows: 3,140
## Columns: 3
## $ mmwr_week <chr> "2019-01", "2019-01", "2019-01", "2019-01", "2019-01", "2019…
## $ category  <chr> "cardiovascular", "digestive", "infectious", "infectious res…
## $ count     <chr> "843", "441", "299", "240", "691", "4090", "427", "220", "52…
colSums(is.na(allcause_state_byweek))
## mmwr_week  category     count 
##         0         0         0

Deficit Function Used

deficit_plot <- function(fit, 
                        title     = "", 
                        ylim      = NULL,
                        show.data = TRUE,
                        alpha     = 0.05){

  requireNamespace("ggplot2")

  z <- qnorm(1 - alpha/2)


  p <- with(fit, data.frame(date = date,
                        y        = (observed - expected)/expected,
                        increase = fitted,
                        sd       = sd,
                        se       = se)) %>%
    ggplot(aes(date, y)) +
    geom_ribbon(aes(ymin = increase - z * se, ymax = increase + z * se), alpha = 0.5, fill = "#3366FF") + #needs to be edited
    geom_hline(yintercept = 0) +
    ylab("% Change from expected count") +
    scale_y_continuous(labels = scales::percent) +
    scale_x_date(date_labels = "%b %Y") +
    ggtitle(title) +
    xlab("Date") 

  if(show.data) p <- p + geom_point(alpha = 0.3)

  if(!is.null(ylim)) p <- p + coord_cartesian(ylim = ylim)

  return(p + geom_line(aes(y = increase), size=0.70, col="#3366FF"))
}
deficit_cumulative <- function(fit, start, end){
  if (!"curve_fit" %in% attr(fit, "type"))
    stop("This is not the correct deficit_model fit, needs curve fit.")

  ind             <- which(fit$date %in% seq(start, end, by = "day"))
  n               <- length(ind) # returns the # of elements in the ind object
  A               <- matrix(1, n, n)
  A[upper.tri(A)] <- 0                                         # It sets all the upper triangular elements of matrix A (excluding the diagonal) to 0
  A               <- sweep(A, 2, fit$expected[ind], FUN = "*") # This line multiplies each column (2 means columns) of the matrix A by a corresponding value from fit$expected[ind]
                                                               # So, now it's scaling each column of A by expected values from the model fit (fit$expected[ind])
  fhat <- matrix(fit$fitted[ind], ncol = 1)                    # creates a column vector (matrix with 1 column) called fhat from the values in fit$fitted[ind]

 
   # Below, this code is constructing a data frame with:
  
  fit_deficit   <- A %*% fhat                                    # multiply design matrix A by fitted values fhat to get cumulative modeled effect.
  obs_deficit   <- cumsum(fit$observed[ind] - fit$expected[ind])
  varcov       <- fit$x[ind,] %*% fit$betacov %*% t(fit$x[ind,]) # Compute variance-covariance matrix of the linear predictor based on model design matrix x and estimated betacov
  diag(varcov) <- fit$se[ind]^2                                  # Replace diagonal entries of varcov with squared standard errors — probably to correct for over/under-dispersion or model misspecification.
  fitted_se    <- sqrt(diag(A %*%  varcov %*% t(A)))
  sd           <- sqrt(diag(A %*% (fit$cov[ind,ind] + diag(fit$log_expected_se[ind]^2)) %*% t(A))) # this computes the variance of the observed cumulative deficit.
                                                                                                   # It adds fit$cov (covariance of observed data) and extra uncertainty from fit$log_expected_se.
  data.frame(date     = fit$date[ind], # The time points
             observed = obs_deficit,   # Cumulative observed - expected difference
             sd       = sd,            # Modeled cumulative deficit/excess from the model.
             fitted   = fit_deficit,   # Standard deviation of the observed cumulative deficit.
             se       = fitted_se)     # Standard error of the fitted cumulative deficit.
}

Deficit Model fits for the different categories

Cardio

exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 2, 20), 
                     by = "day")

# cardiovascular data
cardio <- all_cause_state_by_week %>%
  filter(category == "cardiovascular") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)
3
## [1] 3
# deficit model with COVID-19 dates excluded
cardio_deficit <- compute_expected(cardio, 
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 8.85.
cardio_deficit
## # A tibble: 312 × 12
##    MMWRyear MMWRweek category       weekly_count_statelevel date      
##       <int>    <int> <chr>          <chr>                   <date>    
##  1     2019        2 cardiovascular 1233                    2019-01-06
##  2     2019        3 cardiovascular 1160                    2019-01-13
##  3     2019        4 cardiovascular 1201                    2019-01-20
##  4     2019        5 cardiovascular 1157                    2019-01-27
##  5     2019        6 cardiovascular 1286                    2019-02-03
##  6     2019        7 cardiovascular 1261                    2019-02-10
##  7     2019        8 cardiovascular 1141                    2019-02-17
##  8     2019        9 cardiovascular 1156                    2019-02-24
##  9     2019       10 cardiovascular 1208                    2019-03-03
## 10     2019       11 cardiovascular 1232                    2019-03-10
## # ℹ 302 more rows
## # ℹ 7 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>,
## #   population <dbl>, outcome <dbl>, log_expected_se <dbl>, expected <dbl>,
## #   excluded <lgl>
cardio_deficit <- cardio_deficit %>%
  mutate(date = as.Date(paste(MMWRyear, MMWRweek, 1, sep = "-"), "%Y-%U-%u"))
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `date = as.Date(paste(MMWRyear, MMWRweek, 1, sep = "-"),
##   "%Y-%U-%u")`.
## Caused by warning in `strptime()`:
## ! (0-based) yday 369 in year 2020 is invalid
combined_plot_1 <- ggplot(cardio_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Cardiovascular Hospitalizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
combined_plot_1
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("combined_plot_1.png", plot = combined_plot_1, width = 8, height = 6, dpi = 300)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
fit_x <- cardio_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = as.Date("2019-01-14"),
               end = as.Date("2024-12-31"),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_x <- deficit_plot(fit_x, 
            title = "Deviations from Expected Counts for Cardio in MA"
)
deficit_plot_x

fit_x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-09-23 2019-12-02       8.144139       8.669496    0.07650853    12062
## 2 2020-01-20 2020-02-03       8.276251       8.703193    0.14678711     3343
## 3 2020-03-02 2021-02-15       7.265876       8.874462    0.03594969    49893
## 4 2021-12-27 2022-01-31       8.304721       9.024090    0.10569035     6709
## 5 2023-05-01 2023-05-01       9.677494       9.573722    0.26665497     1303
##    expected     deficit        sd
## 1 12840.088   -778.0880 113.31411
## 2  3515.453   -172.4533  59.29126
## 3 60938.774 -11045.7736 246.85780
## 4  7290.145   -581.1446  85.38234
## 5  1289.028     13.9721  35.90303
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_1 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Cardio",
       title    = "Cumulative Deficit Hospitilization for Cardio Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_1

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_1_adj <- combined_plot_1 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_x_adj <- deficit_plot_x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_1_adj <- cumulative_deficit_plot_1 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_1 <- combined_plot_1_adj /
                       deficit_plot_x_adj /
                       cumulative_deficit_plot_1_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_1
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("Cardio_combined_figure.png",
       plot = final_combined_plot_1,
       width = 11,      
       height = 10.5,    
       dpi = 300,
       bg = "white")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("cardio_combined_figure.png",
       plot = final_combined_plot_1,
       width = 11,       # good for standard document width
       height = 10.5,    # leaves space for three stacked plots
       dpi = 300,
       bg = "white")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

Digestive

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 2, 27), 
                     by = "day")

# digestive data
digestive <- all_cause_state_by_week %>%
  filter(category == "digestive") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
digestive_deficit <- compute_expected(digestive, 
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 4.33.
digestive_deficit
## # A tibble: 312 × 12
##    MMWRyear MMWRweek category  weekly_count_statelevel date      
##       <int>    <int> <chr>     <chr>                   <date>    
##  1     2019        2 digestive 535                     2019-01-06
##  2     2019        3 digestive 570                     2019-01-13
##  3     2019        4 digestive 553                     2019-01-20
##  4     2019        5 digestive 576                     2019-01-27
##  5     2019        6 digestive 609                     2019-02-03
##  6     2019        7 digestive 599                     2019-02-10
##  7     2019        8 digestive 602                     2019-02-17
##  8     2019        9 digestive 562                     2019-02-24
##  9     2019       10 digestive 550                     2019-03-03
## 10     2019       11 digestive 622                     2019-03-10
## # ℹ 302 more rows
## # ℹ 7 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>,
## #   population <dbl>, outcome <dbl>, log_expected_se <dbl>, expected <dbl>,
## #   excluded <lgl>
# Plotting
combined_plot_2 <- ggplot(digestive_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly digestive ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()

combined_plot_2

fit_2x <- digestive_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_2x <- deficit_plot(fit_2x, 
            title = "Deviations from Expected Counts for Digestive in MA"
)
deficit_plot_2x

fit_2x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-09-15 2019-10-13       3.986860       4.261767    0.07956446     2684
## 2 2020-03-01 2020-09-06       3.519113       4.442692    0.03432838    13267
## 3 2020-10-11 2021-03-14       3.714835       4.193817    0.03680019    11504
## 4 2021-12-05 2022-01-30       3.759757       4.186787    0.05877984     4556
## 5 2022-05-01 2022-05-22       4.368984       4.646495    0.09288423     2353
## 6 2022-09-18 2022-10-02       4.263148       4.543836    0.10606203     1722
##    expected    deficit        sd
## 1  2869.071  -185.0706  53.56371
## 2 16748.877 -3481.8769 129.41745
## 3 12987.299 -1483.2986 113.96183
## 4  5073.467  -517.4671  71.22827
## 5  2502.459  -149.4588  50.02458
## 6  1835.377  -113.3774  42.84131
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_2x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_2 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Digestive",
       title    = "Cumulative Deficit Hospitilization for Digestive Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_2

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_2_adj <- combined_plot_2 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_2x_adj <- deficit_plot_2x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_2_adj <- cumulative_deficit_plot_2 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_2 <- combined_plot_2_adj /
                       deficit_plot_2x_adj /
                       cumulative_deficit_plot_2_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_2

ggsave("Digestive_combined_figure.png",
       plot = final_combined_plot_2,
       width = 11,      
       height = 10.5,    
       dpi = 300,
       bg = "white")

Infectious

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 6, 5), 
                     by = "day")

# infectious data
infectious <- all_cause_state_by_week %>%
  filter(category == "infectious") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# deficit model with COVID-19 dates excluded
infectious_deficit <- compute_expected(infectious,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 2.66.
# Plotting
combined_plot_3 <- ggplot(infectious_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Infectious ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_3

fit_3x <- infectious_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

deficit_plot_3x <- deficit_plot(fit_3x, 
            title = "Deviations from Expected Counts for Infectious in MA"
)
deficit_plot_3x

fit_3x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-05-03 2020-08-02       2.384625       2.718529    0.03797626     4495
## 2 2020-10-04 2021-05-30       2.231734       2.726452    0.02405327    10517
## 3 2022-01-23 2022-02-27       2.530161       2.785967    0.05872481     2044
## 4 2022-12-25 2023-01-22       2.565316       2.820316    0.06472516     1727
## 5 2023-06-11 2023-07-02       2.519639       2.766860    0.07167585     1357
## 6 2024-10-20 2024-12-15       2.525210       2.768028    0.04779399     3060
##    expected    deficit        sd
## 1  5124.405  -629.4047  71.58495
## 2 12848.349 -2331.3491 113.35056
## 3  2250.654  -206.6542  47.44106
## 4  1898.669  -171.6689  43.57372
## 5  1490.146  -133.1456  38.60240
## 6  3354.243  -294.2430  57.91583
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_3x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_3 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Infectious",
       title    = "Cumulative Deficit Hospitilization for Infectious Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_3

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_3_adj <- combined_plot_3 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_3x_adj <- deficit_plot_3x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_3_adj <- cumulative_deficit_plot_3 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_3 <- combined_plot_3_adj /
                       deficit_plot_3x_adj /
                       cumulative_deficit_plot_3_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_3

ggsave("Infectious_combined_figure.png",
       plot = final_combined_plot_3,
       width = 11,       
       height = 10.5,    
       dpi = 300,
       bg = "white")

infectious_respiratory

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 6, 12), 
                     by = "day")

# infectious respiratory data
infectious_respiratory <- all_cause_state_by_week %>%
  filter(category == "infectious respiratory") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
infectious_respiratory_deficit <- compute_expected(infectious_respiratory,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 0.699.
# Plotting
combined_plot_4 <- ggplot(infectious_respiratory_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Infectious Respiratory ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_4

# -- Fitting excess model
fit_4x <- infectious_respiratory_deficit %>% 
  deficit_model(exclude = exclude_dates,
                start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_4x <- deficit_plot(fit_4x, 
            title = "Deviations from Expected Counts for Infectious Respiratory in MA"
)
deficit_plot_4x

# -- Intervals of inordinate mortality found by the excess model
fit_4x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-10-13 2019-12-01      0.6340876      1.2580696    0.03417564      683
## 2 2020-03-29 2021-05-30      0.1241042      0.8568669    0.01013141     1036
## 3 2021-10-31 2022-03-20      0.4714432      1.3271874    0.02166535     1333
## 4 2023-01-15 2023-02-26      0.6535837      1.2568241    0.03651720      616
## 5 2023-04-09 2023-04-23      0.4208683      0.7014588    0.04167251      170
## 6 2024-11-03 2024-11-10      0.6424430      1.0037101    0.06105179      173
##    expected     deficit       sd
## 1 1355.1149  -672.11494 36.81189
## 2 7152.9723 -6116.97228 84.57525
## 3 3752.6065 -2419.60646 61.25852
## 4 1184.5517  -568.55175 34.41732
## 5  283.3380  -113.33805 16.83265
## 6  270.2837   -97.28366 16.44031
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_4x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

cumulative_deficit_plot_4 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-infectious respiratory",
       title    = "Cumulative Deficit Hospitilization for infectious respiratory Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_4

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_4_adj <- combined_plot_4 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_4x_adj <- deficit_plot_4x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_4_adj <- cumulative_deficit_plot_4 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_4 <- combined_plot_4_adj /
                       deficit_plot_4x_adj /
                       cumulative_deficit_plot_4_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_4

ggsave("Infectious Respiratory_combined_figure.png",
       plot = final_combined_plot_4,
       width = 11,       # good for standard document width
       height = 10.5,    # leaves space for three stacked plots
       dpi = 300,
       bg = "white")

injury

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 4, 10), 
                     by = "day")

# injury data
injury <- all_cause_state_by_week %>%
  filter(category == "injury") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
injury_deficit <- compute_expected(injury,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 7.86.
# Plotting
combined_plot_5 <- ggplot(injury_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Injury ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_5

# -- Fitting excess model
fit_5x <- injury_deficit %>%                                                     
  deficit_model(exclude = exclude_dates,
                start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_5x <- excess_plot(fit_5x, 
            title = "Deviations from Expected Counts for Injuries in MA"
)
deficit_plot_5x

# -- Intervals of inordinate mortality found by the excess model
fit_5x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-10-20 2019-11-24       6.756174       7.279676    0.09492698     5458
## 2 2020-01-12 2020-09-13       6.349953       7.623419    0.03965819    30779
## 3 2020-10-04 2021-02-07       6.578836       7.492868    0.05411984    16830
## 4 2021-03-07 2021-04-11       7.120101       7.580370    0.09686766     5752
## 5 2021-05-02 2021-05-16       7.600386       7.944070    0.14023944     3070
## 6 2021-07-25 2021-08-08       7.979167       8.372155    0.14396843     3223
## 7 2021-12-19 2022-01-16       7.140402       7.682093    0.10682282     4807
## 8 2022-07-31 2022-08-21       8.268205       8.668967    0.12687117     4453
##    expected    deficit        sd
## 1  5880.914  -422.9136  76.68712
## 2 36951.644 -6172.6438 192.22810
## 3 19168.280 -2338.2799 138.44956
## 4  6123.830  -371.8301  78.25490
## 5  3208.823  -138.8232  56.64648
## 6  3381.738  -158.7384  58.15272
## 7  5171.673  -364.6727  71.91434
## 8  4668.838  -215.8381  68.32890
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_5x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_5 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Injury",
       title    = "Cumulative Deficit Hospitilization for Injury Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_5

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_5_adj <- combined_plot_5 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_5x_adj <- deficit_plot_5x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_5_adj <- cumulative_deficit_plot_5 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_5 <- combined_plot_5_adj /
                       deficit_plot_5x_adj /
                       cumulative_deficit_plot_5_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_5

ggsave("Injury_combined_figure.png",
       plot = final_combined_plot_5,
       width = 11,       # good for standard document width
       height = 10.5,    # leaves space for three stacked plots
       dpi = 300,
       bg = "white")

pneumonia

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2022, 5, 7), 
                     by = "day")

# pneumonia data
pneumonia <- all_cause_state_by_week %>%
  filter(category == "pneumonia") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
pneumonia_deficit <- compute_expected(pneumonia,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 2.32.
# Plotting
combined_plot_6 <- ggplot(pneumonia_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Pneumonia ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_6

# -- Fitting excess model
fit_6x <- pneumonia_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_6x <- deficit_plot(fit_6x, 
            title = "Deviations from Expected Counts for Pneumonia in MA"
)
deficit_plot_6x

# -- Intervals of inordinate mortality found by the excess model
fit_6x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-10-27 2019-12-01       2.532637       2.967751    0.06061043     2046
## 2 2020-04-26 2021-07-04       1.489661       2.618500    0.01756974    12636
## 3 2021-09-05 2022-05-22       1.947460       2.813849    0.02345135     9964
## 4 2022-07-17 2022-08-28       1.694437       2.049678    0.04663405     1597
## 5 2022-11-06 2022-11-13       2.536350       2.922000    0.10416801      683
## 6 2023-01-22 2023-02-19       2.407862       2.978610    0.06651675     1621
## 7 2023-04-23 2023-05-21       2.294970       2.662542    0.06288869     1545
## 8 2024-01-28 2024-02-11       2.619286       2.974941    0.08581986     1058
##     expected    deficit        sd
## 1  2397.5084  -351.5084  48.96436
## 2 22211.3287 -9575.3287 149.03466
## 3 14396.7975 -4432.7975 119.98666
## 4  1931.8138  -334.8138  43.95240
## 5   786.8496  -103.8496  28.05084
## 6  2005.2341  -384.2341  44.77984
## 7  1792.4540  -247.4540  42.33738
## 8  1201.6586  -143.6586  34.66495
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_6x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_6 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-pneumonia",
       title    = "Cumulative Deficit Hospitilization for pneumonia Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_6

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_6_adj <- combined_plot_6 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_6x_adj <- deficit_plot_6x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_6_adj <- cumulative_deficit_plot_6 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_6 <- combined_plot_6_adj /
                       deficit_plot_6x_adj /
                       cumulative_deficit_plot_6_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_6

ggsave("Pneumonia_combined_figure.png",
       plot = final_combined_plot_6,
       width = 11,       
       height = 10.5,    
       dpi = 300,
       bg = "white")

Renal

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 2, 27), 
                     by = "day")

# renal data
renal <- all_cause_state_by_week %>%
  filter(category == "renal") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
renal_deficit <- compute_expected(renal,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 2.42.
# Plotting
combined_plot_7 <- ggplot(renal_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Renal ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_7

# -- Fitting excess model
fit_7x <- renal_deficit %>% 
  deficit_model(exclude = exclude_dates,
                start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_7x <- deficit_plot(fit_7x, 
            title = "Deviations from Expected Counts for Renal in MA"
)
deficit_plot_7x

# -- Intervals of inordinate mortality found by the excess model
fit_7x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-11-10 2019-12-01       1.914332       2.068933    0.06198012     1031
## 2 2020-03-01 2020-07-12       1.798840       2.309561    0.02928593     4844
## 3 2020-09-06 2020-10-04       2.170195       2.326802    0.05879006     1461
## 4 2020-11-01 2021-02-28       1.966115       2.218355    0.03025440     4765
## 5 2021-12-26 2022-01-02       2.213272       2.330685    0.09303280      596
## 6 2022-05-15 2022-05-29       2.364289       2.607537    0.08034593      955
##    expected     deficit       sd
## 1 1114.2633   -83.26328 33.38058
## 2 6219.2915 -1375.29147 78.86248
## 3 1566.4295  -105.42948 39.57814
## 4 5376.3201  -611.32014 73.32339
## 5  627.6176   -31.61763 25.05230
## 6 1053.2544   -98.25443 32.45388
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_7x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_7 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-renal",
       title    = "Cumulative Deficit Hospitilization for renal Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_7

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_7_adj <- combined_plot_7 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_7x_adj <- deficit_plot_7x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_7_adj <- cumulative_deficit_plot_7 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_7 <- combined_plot_7_adj /
                       deficit_plot_7x_adj /
                       cumulative_deficit_plot_7_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_7

ggsave("Renal_combined_figure.png",
       plot = final_combined_plot_7,
       width = 11,       # good for standard document width
       height = 10.5,    # leaves space for three stacked plots
       dpi = 300,
       bg = "white")

Respiratory

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 6, 12), 
                     by = "day")

# respiratory data
respiratory <- all_cause_state_by_week %>%
  filter(category == "respiratory") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
respiratory_deficit <- compute_expected(respiratory,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 3.89.
# Plotting
combined_plot_8 <- ggplot(respiratory_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Respiratory ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_8

# -- Fitting excess model
fit_8x <- respiratory_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_8x <- deficit_plot(fit_8x, 
            title = "Deviations from Expected Counts for Respiratory in MA"
)
deficit_plot_8x

# -- Intervals of inordinate mortality found by the excess model
fit_8x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-07-21 2019-07-21       2.978262       3.366797    0.15813126      401
## 2 2019-10-27 2019-11-24       3.920017       4.324783    0.08015053     2639
## 3 2020-03-08 2021-05-30       2.675694       4.168126    0.02182343    23417
## 4 2021-11-21 2022-03-27       3.632627       4.441859    0.04166917     9293
## 5 2023-04-23 2023-05-07       4.186401       4.696742    0.10783183     1691
##     expected      deficit        sd
## 1   453.3133    -52.31326  21.29115
## 2  2911.4933   -272.49331  53.95826
## 3 36478.3943 -13061.39427 190.99318
## 4 11363.1786  -2070.17861 106.59821
## 5  1897.1403   -206.14030  43.55617
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_8x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_8 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-respiratory",
       title    = "Cumulative Deficit Hospitilization for respiratory Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_8

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_8_adj <- combined_plot_8 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_8x_adj <- deficit_plot_8x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_8_adj <- cumulative_deficit_plot_8 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_8 <- combined_plot_8_adj /
                       deficit_plot_8x_adj /
                       cumulative_deficit_plot_8_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_8

ggsave("Respiratory_combined_figure.png",
       plot = final_combined_plot_8,
       width = 11,      
       height = 10.5,   
       dpi = 300,
       bg = "white")

other

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 6, 5), 
                     by = "day")

# other data
other <- all_cause_state_by_week %>%
  filter(category == "other") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
other_deficit <- compute_expected(other,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 41.3.
# Plotting
combined_plot_9 <- ggplot(other_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Other ED Hospitilizations in MA",
       x = "Date",
       y = "Weekly Hospitilizations",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_9

# -- Fitting excess model
fit_9x <-  other_deficit %>% 
  deficit_model(exclude = exclude_dates,
                start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_9x <- deficit_plot(fit_9x, 
            title = "Deviations from Expected Counts for Other in MA"
)
deficit_plot_9x

# -- Intervals of inordinate mortality found by the excess model
fit_9x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-03-01 2021-05-23       35.68635       41.60681    0.06895008   312318
## 2 2021-08-22 2021-10-10       40.63174       42.79761    0.19933061    43766
## 3 2021-11-21 2022-02-13       36.00309       40.26895    0.15167805    63018
## 4 2022-03-06 2022-04-03       41.10001       42.97170    0.25264781    27669
## 5 2022-04-24 2022-05-08       42.59682       44.33594    0.33130392    17206
##    expected     deficit       sd
## 1 364132.31 -51814.3141 603.4338
## 2  46098.94  -2332.9402 214.7066
## 3  70484.75  -7466.7499 265.4896
## 4  28929.04  -1260.0429 170.0854
## 5  17908.48   -702.4761 133.8226
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_9x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_9 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-other",
       title    = "Cumulative Deficit Hospitilization for other Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_9

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_9_adj <- combined_plot_9 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_9x_adj <- deficit_plot_9x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_9_adj <- cumulative_deficit_plot_9 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_9 <- combined_plot_9_adj /
                       deficit_plot_9x_adj /
                       cumulative_deficit_plot_9_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_9

ggsave("Other_combined_figure.png",
       plot = final_combined_plot_9,
       width = 11,       
       height = 10.5,    
       dpi = 300,
       bg = "white")

Different Approuch

library(ggplot2)
library(patchwork)
library(cowplot)
library(scales)

# Common scale settings
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(x_limits[1], x_limits[2], by = "1 year")
category_labels <- c("Cardiovascular", "Digestive", "Infectious", "Infectious Respiratory",
                     "Injury", "Pneumonia", "Renal", "Respiratory", "Other")

# ----- CLEANERS -----
clean_combined <- function(plots, labels) {
  mapply(function(p, label) {
    p +
      labs(title = label, x = NULL) +
      scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
      theme_minimal() +
      theme(legend.position = "none",
            plot.title = element_text(hjust = 0.5, size = 11),
            axis.text.x = element_text(size = 8))
  }, plots, labels, SIMPLIFY = FALSE)
}

clean_deficit <- function(plot_list, labels) {
  mapply(function(p, label) {
    p +
      labs(title = label, x = NULL) +
      scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
      theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5, size = 11),
            axis.text.x = element_text(size = 8))
  }, plot_list, labels, SIMPLIFY = FALSE)
}

clean_cumulative <- function(plots, labels) {
  mapply(function(p, label) {
    p +
      labs(title = label, x = NULL, y = NULL, subtitle = NULL) +
      scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
      scale_y_continuous(labels = comma) +
      theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5, size = 11),
            axis.text = element_text(size = 8))
  }, plots, labels, SIMPLIFY = FALSE)
}

# ----- APPLY CLEANING -----
combined_clean <- clean_combined(list(
  combined_plot_1, combined_plot_2, combined_plot_3,
  combined_plot_4, combined_plot_5, combined_plot_6,
  combined_plot_7, combined_plot_8, combined_plot_9
), category_labels)

deficit_clean <- clean_deficit(list(
  deficit_plot_x, deficit_plot_2x, deficit_plot_3x,
  deficit_plot_4x, deficit_plot_5x, deficit_plot_6x,
  deficit_plot_7x, deficit_plot_8x, deficit_plot_9x
), category_labels)
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_clean <- clean_cumulative(list(
  cumulative_deficit_plot_1, cumulative_deficit_plot_2, cumulative_deficit_plot_3,
  cumulative_deficit_plot_4, cumulative_deficit_plot_5, cumulative_deficit_plot_6,
  cumulative_deficit_plot_7, cumulative_deficit_plot_8, cumulative_deficit_plot_9
), category_labels)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
# ----- LEGEND FOR COMBINED -----
legend_shared <- get_legend(
  combined_plot_1 +
    theme(legend.position = "bottom",
          legend.title = element_text(size = 10),
          legend.text = element_text(size = 9))
)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.
# ----- FINAL COMBINED PLOT (balanced, readable) -----
combined_final <- wrap_plots(combined_clean, ncol = 3) +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

combined_final <- combined_final +
  plot_annotation(
    title = "Weekly Observed vs Expected Hospitalizations by Cause",
    theme = theme(plot.title = element_text(hjust = 0.5, size = 15))
  )

ggsave("combined_all_categories.png",
       plot = combined_final,
       width = 12, height = 10, dpi = 300, bg = "white")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
deficit_final <- wrap_plots(deficit_clean, ncol = 3) +
  plot_annotation(
    title = "Deviations from Expected Weekly Admissions",
    theme = theme(plot.title = element_text(hjust = 0.5, size = 14))
  )

ggsave("deficit_all_categories.png",
       plot = deficit_final,
       width = 11, height = 9.5, dpi = 300, bg = "white")


cumulative_grid <- wrap_plots(cumulative_clean, ncol = 3)

y_label <- ggdraw() +
  draw_label("Cumulative Deficit in Hospitalizations", angle = 90, size = 12)

cumulative_final <- plot_grid(
  y_label, cumulative_grid,
  rel_widths = c(0.04, 0.96),
  nrow = 1
)

cumulative_final <- cumulative_final +
  plot_annotation(
    title = "Cumulative Deficit in Hospital Admissions by Cause (2020–2024)",
    theme = theme(plot.title = element_text(hjust = 0.5, size = 15))
  )

ggsave("cumulative_all_categories.png",
       plot = cumulative_final,
       width = 12, height = 10, dpi = 300, bg = "white")

combined_final
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

deficit_final

cumulative_final

ggsave("deficit_plot_x.png", plot = deficit_plot_x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_2x.png", plot = deficit_plot_2x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_3x.png", plot = deficit_plot_3x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_4x.png", plot = deficit_plot_4x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_5x.png", plot = deficit_plot_5x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_6x.png", plot = deficit_plot_6x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_7x.png", plot = deficit_plot_7x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_8x.png", plot = deficit_plot_8x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_9x.png", plot = deficit_plot_9x, width = 8, height = 6, dpi = 300)

Cumulative deficit interactive plots

cumulative_deficit <- deficit_cumulative(fit_x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for Cardio Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-Cardio", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_2x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_2 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Digestive",
       title    = "Cumulative Deficit Hospitilization for Digestive Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_2

cumulative_deficit <- deficit_cumulative(fit_2x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for Digestive Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-Digestive", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#3

cumulative_deficit <- deficit_cumulative(fit_3x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for Infectious Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-Infectious", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#4

cumulative_deficit <- deficit_cumulative(fit_4x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for infectious respiratory Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-infectious respiratory", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#5

cumulative_deficit <- deficit_cumulative(fit_5x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for injury Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-injury", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#6

cumulative_deficit <- deficit_cumulative(fit_6x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for pneumonia Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-pneumonia", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#7

cumulative_deficit <- deficit_cumulative(fit_7x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for renal Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-renal", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#8

cumulative_deficit <- deficit_cumulative(fit_8x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for respiratory Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-respiratory", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#9

cumulative_deficit <- deficit_cumulative(fit_9x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for other Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-other", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly
ggsave("cumulative_deficit_plot_1.png", plot = cumulative_deficit_plot_1, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_2.png", plot = cumulative_deficit_plot_2, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_3.png", plot = cumulative_deficit_plot_3, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_4.png", plot = cumulative_deficit_plot_4, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_5.png", plot = cumulative_deficit_plot_5, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_6.png", plot = cumulative_deficit_plot_6, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_7.png", plot = cumulative_deficit_plot_7, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_8.png", plot = cumulative_deficit_plot_8, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_9.png", plot = cumulative_deficit_plot_9, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_1.png", plot = cumulative_deficit_plot_1, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_2.png", plot = cumulative_deficit_plot_2, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_3.png", plot = cumulative_deficit_plot_3, width = 5, height = 3, dpi = 300)

ggsave("cumulative_deficit_plot_4.png", plot = cumulative_deficit_plot_4, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_5.png", plot = cumulative_deficit_plot_5, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_6.png", plot = cumulative_deficit_plot_6, width = 5, height = 3, dpi = 300)

ggsave("cumulative_deficit_plot_7.png", plot = cumulative_deficit_plot_7, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_8.png", plot = cumulative_deficit_plot_8, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_9.png", plot = cumulative_deficit_plot_9, width = 5, height = 3, dpi = 300)
# new 'category' column added
cardiovascular <- fit_x$detected_intervals %>%
  mutate(category = "cardiovascular")
digestive <- fit_2x$detected_intervals %>%
  mutate(category = "digestive")
infectious <- fit_3x$detected_intervals %>%
  mutate(category = "infectious")
infectious_respiratory <- fit_4x$detected_intervals %>%
  mutate(category = "infectious respiratory")
injury <- fit_5x$detected_intervals %>%
  mutate(category = "injury")
pneumonia <- fit_6x$detected_intervals %>%
  mutate(category = "pneumonia")
renal <- fit_7x$detected_intervals %>%
  mutate(category = "renal")
respiratory <- fit_8x$detected_intervals %>%
  mutate(category = "respiratory")
other <- fit_9x$detected_intervals %>%
  mutate(category = "other")

combined_intervals <- bind_rows(
  cardiovascular,
  digestive,
  infectious,
  infectious_respiratory,
  injury,
  pneumonia,
  renal,
  respiratory,
  other
)

combined_intervals <- combined_intervals %>%
  select(category, everything())

combined_intervals <- combined_intervals %>%
  mutate(across(c(where(is.numeric)), ~round(., 1)))

kable_table <- combined_intervals %>%
  kable(
    format = "html",
    caption = "Detected Intervals by Disease Category",
    align = c("l", rep("c", 9)),
    col.names = c(
      "Category", "Start Date", "End Date", "Observed Count Rate", 
      "Expected Count Rate", "SD Count Rate", "Observed Counts", 
      "Expected Counts", "Deficit", "SD"
    ),
    row.names = FALSE
  )

kable_table
Detected Intervals by Disease Category
Category Start Date End Date Observed Count Rate Expected Count Rate SD Count Rate Observed Counts Expected Counts Deficit SD
cardiovascular 2019-09-23 2019-12-02 8.1 8.7 0.1 12062 12840.1 -778.1 113.3
cardiovascular 2020-01-20 2020-02-03 8.3 8.7 0.1 3343 3515.5 -172.5 59.3
cardiovascular 2020-03-02 2021-02-15 7.3 8.9 0.0 49893 60938.8 -11045.8 246.9
cardiovascular 2021-12-27 2022-01-31 8.3 9.0 0.1 6709 7290.1 -581.1 85.4
cardiovascular 2023-05-01 2023-05-01 9.7 9.6 0.3 1303 1289.0 14.0 35.9
digestive 2019-09-15 2019-10-13 4.0 4.3 0.1 2684 2869.1 -185.1 53.6
digestive 2020-03-01 2020-09-06 3.5 4.4 0.0 13267 16748.9 -3481.9 129.4
digestive 2020-10-11 2021-03-14 3.7 4.2 0.0 11504 12987.3 -1483.3 114.0
digestive 2021-12-05 2022-01-30 3.8 4.2 0.1 4556 5073.5 -517.5 71.2
digestive 2022-05-01 2022-05-22 4.4 4.6 0.1 2353 2502.5 -149.5 50.0
digestive 2022-09-18 2022-10-02 4.3 4.5 0.1 1722 1835.4 -113.4 42.8
infectious 2020-05-03 2020-08-02 2.4 2.7 0.0 4495 5124.4 -629.4 71.6
infectious 2020-10-04 2021-05-30 2.2 2.7 0.0 10517 12848.3 -2331.3 113.4
infectious 2022-01-23 2022-02-27 2.5 2.8 0.1 2044 2250.7 -206.7 47.4
infectious 2022-12-25 2023-01-22 2.6 2.8 0.1 1727 1898.7 -171.7 43.6
infectious 2023-06-11 2023-07-02 2.5 2.8 0.1 1357 1490.1 -133.1 38.6
infectious 2024-10-20 2024-12-15 2.5 2.8 0.0 3060 3354.2 -294.2 57.9
infectious respiratory 2019-10-13 2019-12-01 0.6 1.3 0.0 683 1355.1 -672.1 36.8
infectious respiratory 2020-03-29 2021-05-30 0.1 0.9 0.0 1036 7153.0 -6117.0 84.6
infectious respiratory 2021-10-31 2022-03-20 0.5 1.3 0.0 1333 3752.6 -2419.6 61.3
infectious respiratory 2023-01-15 2023-02-26 0.7 1.3 0.0 616 1184.6 -568.6 34.4
infectious respiratory 2023-04-09 2023-04-23 0.4 0.7 0.0 170 283.3 -113.3 16.8
infectious respiratory 2024-11-03 2024-11-10 0.6 1.0 0.1 173 270.3 -97.3 16.4
injury 2019-10-20 2019-11-24 6.8 7.3 0.1 5458 5880.9 -422.9 76.7
injury 2020-01-12 2020-09-13 6.3 7.6 0.0 30779 36951.6 -6172.6 192.2
injury 2020-10-04 2021-02-07 6.6 7.5 0.1 16830 19168.3 -2338.3 138.4
injury 2021-03-07 2021-04-11 7.1 7.6 0.1 5752 6123.8 -371.8 78.3
injury 2021-05-02 2021-05-16 7.6 7.9 0.1 3070 3208.8 -138.8 56.6
injury 2021-07-25 2021-08-08 8.0 8.4 0.1 3223 3381.7 -158.7 58.2
injury 2021-12-19 2022-01-16 7.1 7.7 0.1 4807 5171.7 -364.7 71.9
injury 2022-07-31 2022-08-21 8.3 8.7 0.1 4453 4668.8 -215.8 68.3
pneumonia 2019-10-27 2019-12-01 2.5 3.0 0.1 2046 2397.5 -351.5 49.0
pneumonia 2020-04-26 2021-07-04 1.5 2.6 0.0 12636 22211.3 -9575.3 149.0
pneumonia 2021-09-05 2022-05-22 1.9 2.8 0.0 9964 14396.8 -4432.8 120.0
pneumonia 2022-07-17 2022-08-28 1.7 2.0 0.0 1597 1931.8 -334.8 44.0
pneumonia 2022-11-06 2022-11-13 2.5 2.9 0.1 683 786.8 -103.8 28.1
pneumonia 2023-01-22 2023-02-19 2.4 3.0 0.1 1621 2005.2 -384.2 44.8
pneumonia 2023-04-23 2023-05-21 2.3 2.7 0.1 1545 1792.5 -247.5 42.3
pneumonia 2024-01-28 2024-02-11 2.6 3.0 0.1 1058 1201.7 -143.7 34.7
renal 2019-11-10 2019-12-01 1.9 2.1 0.1 1031 1114.3 -83.3 33.4
renal 2020-03-01 2020-07-12 1.8 2.3 0.0 4844 6219.3 -1375.3 78.9
renal 2020-09-06 2020-10-04 2.2 2.3 0.1 1461 1566.4 -105.4 39.6
renal 2020-11-01 2021-02-28 2.0 2.2 0.0 4765 5376.3 -611.3 73.3
renal 2021-12-26 2022-01-02 2.2 2.3 0.1 596 627.6 -31.6 25.1
renal 2022-05-15 2022-05-29 2.4 2.6 0.1 955 1053.3 -98.3 32.5
respiratory 2019-07-21 2019-07-21 3.0 3.4 0.2 401 453.3 -52.3 21.3
respiratory 2019-10-27 2019-11-24 3.9 4.3 0.1 2639 2911.5 -272.5 54.0
respiratory 2020-03-08 2021-05-30 2.7 4.2 0.0 23417 36478.4 -13061.4 191.0
respiratory 2021-11-21 2022-03-27 3.6 4.4 0.0 9293 11363.2 -2070.2 106.6
respiratory 2023-04-23 2023-05-07 4.2 4.7 0.1 1691 1897.1 -206.1 43.6
other 2020-03-01 2021-05-23 35.7 41.6 0.1 312318 364132.3 -51814.3 603.4
other 2021-08-22 2021-10-10 40.6 42.8 0.2 43766 46098.9 -2332.9 214.7
other 2021-11-21 2022-02-13 36.0 40.3 0.2 63018 70484.7 -7466.7 265.5
other 2022-03-06 2022-04-03 41.1 43.0 0.3 27669 28929.0 -1260.0 170.1
other 2022-04-24 2022-05-08 42.6 44.3 0.3 17206 17908.5 -702.5 133.8
library(kableExtra)
html_table <- combined_intervals %>%
  kable("html",
        escape = FALSE,
        caption = "Detected Intervals by Disease Category",
        align = c("l", rep("c", ncol(combined_intervals) - 1)),
        col.names = c(
          "Category", "Start Date", "End Date", "Observed Count Rate", 
          "Expected Count Rate", "SD Count Rate", "Observed Counts", 
          "Expected Counts", "Deficit", "SD"
        )) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                position = "center")

html_table
Detected Intervals by Disease Category
Category Start Date End Date Observed Count Rate Expected Count Rate SD Count Rate Observed Counts Expected Counts Deficit SD
1…1 cardiovascular 2019-09-23 2019-12-02 8.1 8.7 0.1 12062 12840.1 -778.1 113.3
2…2 cardiovascular 2020-01-20 2020-02-03 8.3 8.7 0.1 3343 3515.5 -172.5 59.3
3…3 cardiovascular 2020-03-02 2021-02-15 7.3 8.9 0.0 49893 60938.8 -11045.8 246.9
4…4 cardiovascular 2021-12-27 2022-01-31 8.3 9.0 0.1 6709 7290.1 -581.1 85.4
5…5 cardiovascular 2023-05-01 2023-05-01 9.7 9.6 0.3 1303 1289.0 14.0 35.9
1…6 digestive 2019-09-15 2019-10-13 4.0 4.3 0.1 2684 2869.1 -185.1 53.6
2…7 digestive 2020-03-01 2020-09-06 3.5 4.4 0.0 13267 16748.9 -3481.9 129.4
3…8 digestive 2020-10-11 2021-03-14 3.7 4.2 0.0 11504 12987.3 -1483.3 114.0
4…9 digestive 2021-12-05 2022-01-30 3.8 4.2 0.1 4556 5073.5 -517.5 71.2
5…10 digestive 2022-05-01 2022-05-22 4.4 4.6 0.1 2353 2502.5 -149.5 50.0
6…11 digestive 2022-09-18 2022-10-02 4.3 4.5 0.1 1722 1835.4 -113.4 42.8
1…12 infectious 2020-05-03 2020-08-02 2.4 2.7 0.0 4495 5124.4 -629.4 71.6
2…13 infectious 2020-10-04 2021-05-30 2.2 2.7 0.0 10517 12848.3 -2331.3 113.4
3…14 infectious 2022-01-23 2022-02-27 2.5 2.8 0.1 2044 2250.7 -206.7 47.4
4…15 infectious 2022-12-25 2023-01-22 2.6 2.8 0.1 1727 1898.7 -171.7 43.6
5…16 infectious 2023-06-11 2023-07-02 2.5 2.8 0.1 1357 1490.1 -133.1 38.6
6…17 infectious 2024-10-20 2024-12-15 2.5 2.8 0.0 3060 3354.2 -294.2 57.9
1…18 infectious respiratory 2019-10-13 2019-12-01 0.6 1.3 0.0 683 1355.1 -672.1 36.8
2…19 infectious respiratory 2020-03-29 2021-05-30 0.1 0.9 0.0 1036 7153.0 -6117.0 84.6
3…20 infectious respiratory 2021-10-31 2022-03-20 0.5 1.3 0.0 1333 3752.6 -2419.6 61.3
4…21 infectious respiratory 2023-01-15 2023-02-26 0.7 1.3 0.0 616 1184.6 -568.6 34.4
5…22 infectious respiratory 2023-04-09 2023-04-23 0.4 0.7 0.0 170 283.3 -113.3 16.8
6…23 infectious respiratory 2024-11-03 2024-11-10 0.6 1.0 0.1 173 270.3 -97.3 16.4
1…24 injury 2019-10-20 2019-11-24 6.8 7.3 0.1 5458 5880.9 -422.9 76.7
2…25 injury 2020-01-12 2020-09-13 6.3 7.6 0.0 30779 36951.6 -6172.6 192.2
3…26 injury 2020-10-04 2021-02-07 6.6 7.5 0.1 16830 19168.3 -2338.3 138.4
4…27 injury 2021-03-07 2021-04-11 7.1 7.6 0.1 5752 6123.8 -371.8 78.3
5…28 injury 2021-05-02 2021-05-16 7.6 7.9 0.1 3070 3208.8 -138.8 56.6
6…29 injury 2021-07-25 2021-08-08 8.0 8.4 0.1 3223 3381.7 -158.7 58.2
7…30 injury 2021-12-19 2022-01-16 7.1 7.7 0.1 4807 5171.7 -364.7 71.9
8…31 injury 2022-07-31 2022-08-21 8.3 8.7 0.1 4453 4668.8 -215.8 68.3
1…32 pneumonia 2019-10-27 2019-12-01 2.5 3.0 0.1 2046 2397.5 -351.5 49.0
2…33 pneumonia 2020-04-26 2021-07-04 1.5 2.6 0.0 12636 22211.3 -9575.3 149.0
3…34 pneumonia 2021-09-05 2022-05-22 1.9 2.8 0.0 9964 14396.8 -4432.8 120.0
4…35 pneumonia 2022-07-17 2022-08-28 1.7 2.0 0.0 1597 1931.8 -334.8 44.0
5…36 pneumonia 2022-11-06 2022-11-13 2.5 2.9 0.1 683 786.8 -103.8 28.1
6…37 pneumonia 2023-01-22 2023-02-19 2.4 3.0 0.1 1621 2005.2 -384.2 44.8
7…38 pneumonia 2023-04-23 2023-05-21 2.3 2.7 0.1 1545 1792.5 -247.5 42.3
8…39 pneumonia 2024-01-28 2024-02-11 2.6 3.0 0.1 1058 1201.7 -143.7 34.7
1…40 renal 2019-11-10 2019-12-01 1.9 2.1 0.1 1031 1114.3 -83.3 33.4
2…41 renal 2020-03-01 2020-07-12 1.8 2.3 0.0 4844 6219.3 -1375.3 78.9
3…42 renal 2020-09-06 2020-10-04 2.2 2.3 0.1 1461 1566.4 -105.4 39.6
4…43 renal 2020-11-01 2021-02-28 2.0 2.2 0.0 4765 5376.3 -611.3 73.3
5…44 renal 2021-12-26 2022-01-02 2.2 2.3 0.1 596 627.6 -31.6 25.1
6…45 renal 2022-05-15 2022-05-29 2.4 2.6 0.1 955 1053.3 -98.3 32.5
1…46 respiratory 2019-07-21 2019-07-21 3.0 3.4 0.2 401 453.3 -52.3 21.3
2…47 respiratory 2019-10-27 2019-11-24 3.9 4.3 0.1 2639 2911.5 -272.5 54.0
3…48 respiratory 2020-03-08 2021-05-30 2.7 4.2 0.0 23417 36478.4 -13061.4 191.0
4…49 respiratory 2021-11-21 2022-03-27 3.6 4.4 0.0 9293 11363.2 -2070.2 106.6
5…50 respiratory 2023-04-23 2023-05-07 4.2 4.7 0.1 1691 1897.1 -206.1 43.6
1…51 other 2020-03-01 2021-05-23 35.7 41.6 0.1 312318 364132.3 -51814.3 603.4
2…52 other 2021-08-22 2021-10-10 40.6 42.8 0.2 43766 46098.9 -2332.9 214.7
3…53 other 2021-11-21 2022-02-13 36.0 40.3 0.2 63018 70484.7 -7466.7 265.5
4…54 other 2022-03-06 2022-04-03 41.1 43.0 0.3 27669 28929.0 -1260.0 170.1
5…55 other 2022-04-24 2022-05-08 42.6 44.3 0.3 17206 17908.5 -702.5 133.8
save_kable(html_table, "combined_intervals_table.html")
webshot("combined_intervals_table.html", "combined_intervals_table.png", vwidth = 1200, vheight = 800)

combined_intervals <- combined_intervals %>%
  mutate(across(where(is.numeric), ~ round(., 1)))

table_1 <- combined_intervals %>%
  slice(1:5) %>%
  kbl(
    format = "html",
    caption = "Detected Intervals by Category (1 of 2)",
    align = c("l", rep("c", 9)),
    col.names = c(
      "Category", "Start Date", "End Date", "Observed Count Rate", 
      "Expected Count Rate", "SD Count Rate", "Observed Counts", 
      "Expected Counts", "Deficit", "SD"
    )
  ) %>%
  kable_styling(full_width = FALSE)
save_kable(table_1, file = "detected_intervals_1.png", zoom = 3)
## save_kable will have the best result with magick installed.
table_2 <- combined_intervals %>%
  slice(6:9) %>%
  kbl(
    format = "html",
    caption = "Detected Intervals by Category (2 of 2)",
    align = c("l", rep("c", 9)),
    col.names = c(
      "Category", "Start Date", "End Date", "Observed Count Rate", 
      "Expected Count Rate", "SD Count Rate", "Observed Counts", 
      "Expected Counts", "Deficit", "SD"
    )
  ) %>%
  kable_styling(full_width = FALSE)
save_kable(table_2, file = "detected_intervals_2.png", zoom = 3)
## save_kable will have the best result with magick installed.
table_2 
Detected Intervals by Category (2 of 2)
Category Start Date End Date Observed Count Rate Expected Count Rate SD Count Rate Observed Counts Expected Counts Deficit SD
digestive 2019-09-15 2019-10-13 4.0 4.3 0.1 2684 2869.1 -185.1 53.6
digestive 2020-03-01 2020-09-06 3.5 4.4 0.0 13267 16748.9 -3481.9 129.4
digestive 2020-10-11 2021-03-14 3.7 4.2 0.0 11504 12987.3 -1483.3 114.0
digestive 2021-12-05 2022-01-30 3.8 4.2 0.1 4556 5073.5 -517.5 71.2

Further Analysis

library(ggplot2)
library(scales)

total_deficits <- data.frame(
  category = c("Cardiovascular", "Digestive", "Infectious", "Infectious Respiratory",
               "Injury", "Pneumonia", "Renal", "Respiratory", "Other"),
  total_deficit = c(
    sum(cumulative_deficit_plot_1$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_2$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_3$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_4$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_5$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_6$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_7$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_8$data$observed, na.rm = TRUE),
    sum(cumulative_deficit_plot_9$data$observed, na.rm = TRUE)
  )
)

bar_plot <- ggplot(total_deficits, aes(x = reorder(category, total_deficit), y = total_deficit)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  scale_y_continuous(labels = comma) +
  labs(title = "Total Cumulative Deficits by Category (2020–2024)",
       x = "ICD-10 Category",
       y = "Total Cumulative Deficit in Hospitalizations") +
  coord_flip() +
  theme_minimal(base_size = 12)
bar_plot

ggsave("Total_Cumulative_Deficit_by_Category.png",
       plot = bar_plot,
       width = 8.5,
       height = 6,
       dpi = 300,
       bg = "white")
cardio_deficit     <- tail(deficit_cumulative(fit_x,   start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
digestive_deficit  <- tail(deficit_cumulative(fit_2x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
infectious_deficit <- tail(deficit_cumulative(fit_3x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
inf_resp_deficit   <- tail(deficit_cumulative(fit_4x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
injury_deficit     <- tail(deficit_cumulative(fit_5x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
pneumonia_deficit  <- tail(deficit_cumulative(fit_6x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
renal_deficit      <- tail(deficit_cumulative(fit_7x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
respiratory_deficit<- tail(deficit_cumulative(fit_8x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed
other_deficit      <- tail(deficit_cumulative(fit_9x,  start = make_date(2020, 3, 1), end = make_date(2024, 12, 28)), 1)$observed

deficit_summary <- tibble(
  Category = c("Cardiovascular", "Digestive", "Infectious", "Infectious Respiratory", 
               "Injury", "Pneumonia", "Renal", "Respiratory", "Other"),
  Cumulative_Deficit = c(cardio_deficit, digestive_deficit, infectious_deficit,
                         inf_resp_deficit, injury_deficit, pneumonia_deficit,
                         renal_deficit, respiratory_deficit, other_deficit)
)

bar_plot2 <- ggplot(deficit_summary, aes(x = reorder(Category, -Cumulative_Deficit), 
                            y = Cumulative_Deficit)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  geom_text(aes(label = scales::comma(round(Cumulative_Deficit))), 
            vjust = -0.5, size = 3) +
  scale_y_continuous(labels = scales::comma) +
  labs(
    title = "Cumulative Deficit in Hospitalizations by Category",
    subtitle = "Massachusetts, March 2020 – December 2024",
    x = "ICD-10 Category",
    y = "Cumulative Deficit (Observed - Expected)"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
bar_plot2

ggsave("Cumulative Deficit in Hospitalizations by Category.png",
       plot = bar_plot2,
       width = 8.5,
       height = 6,
       dpi = 300,
       bg = "white")

ICD-10 codes table

library(tidyverse)
library(readxl)
library(gt)
library(stringr)

# Load and clean data
icd10categories <- read_excel("icd10categories_updated_321.xlsx") %>%
  mutate(icd10code = toupper(trimws(icd10code)))

# Function to generate code ranges and extract sort keys
get_code_range_with_sort <- function(codes) {
  codes <- sort(unique(codes))
  letter <- substr(codes, 1, 1)
  num <- as.integer(str_extract(substr(codes, 2, 4), "\\d+"))
  df <- tibble(code = codes, letter = letter, num = num)
  
  ranges <- df %>%
    group_by(letter) %>%
    summarize(
      min_num = min(num),
      max_num = max(num),
      .groups = "drop"
    ) %>%
    mutate(
      range = ifelse(min_num == max_num,
                     paste0(letter, sprintf("%02d", min_num)),
                     paste0(letter, sprintf("%02d", min_num), "–", letter, sprintf("%02d", max_num))),
      sort_key = min_num + (utf8ToInt(letter) * 1000) # sort letter first, then number
    )
  
  list(
    range = paste(ranges$range, collapse = ", "),
    sort_key = min(ranges$sort_key)
  )
}

# Apply function and unpack results
icd10_summary <- icd10categories %>%
  group_by(category_0321update) %>%
  summarize(
    `Number of ICD-10 codes` = n_distinct(icd10code),
    .groups = "drop",
    temp = list(get_code_range_with_sort(icd10code))
  ) %>%
  mutate(
    `ICD-10 code groups` = map_chr(temp, "range"),
    sort_key = map_dbl(temp, "sort_key")
  ) %>%
  select(-temp) %>%
  rename(category = category_0321update) %>%
  arrange(sort_key)  # Sort by ICD-10 code order
## Warning: There were 3 warnings in `summarize()`.
## The first warning was:
## ℹ In argument: `sort_key = min_num + (utf8ToInt(letter) * 1000)`.
## ℹ In group 0: .
## Caused by warning:
## ! There was 1 warning in `mutate()`.
## ℹ In argument: `sort_key = min_num + (utf8ToInt(letter) * 1000)`.
## Caused by warning in `utf8ToInt()`:
## ! argument should be a character vector of length 1
## all but the first element will be ignored
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 2 remaining warnings.
# Create GT table
gt_table3 <- icd10_summary %>%
  select(-sort_key) %>%
  gt() %>%
  tab_header(
    title = md("**Table. ICD-10 Code Groupings and Counts by Category**")
  ) %>%
  cols_width(
    `ICD-10 code groups` ~ px(400)
  ) %>%
  tab_options(table.width = pct(100))
gt_table3
Table. ICD-10 Code Groupings and Counts by Category
category Number of ICD-10 codes ICD-10 code groups
infectious 199 A00–A99, B01–B99
cardiovascular 100 I00–I99
infectious respiratory 13 J00–J22
pneumonia 7 J12–J18
respiratory 78 J19–J99
digestive 25 K20–K92
injury 342 M48–M97, R45, S00–S99, T00–T88, V00–V98, W00–W88, X00–X99, Y00–Y99
renal 29 N00–N28
other 75 O00–O99, Z31–Z33
covid 1 U07
# Save the GT table
gtsave(gt_table3, filename = "icd10_category_groupings_continuous.png")
library(tidyverse)
library(readxl)
library(gt)
library(stringr)

# Load and clean data
icd10categories <- read_excel("icd10categories_updated_321.xlsx") %>%
  mutate(icd10code = toupper(trimws(icd10code)))

# Modified condensation function with relaxed grouping logic
condense_icd10_codes <- function(codes, max_gap = 2) {
  codes <- sort(unique(codes))
  letter_part <- substr(codes, 1, 1)
  num_part <- as.integer(str_extract(substr(codes, 2, 4), "\\d+"))
  
  df <- tibble(code = codes, letter = letter_part, num = num_part) %>%
    arrange(letter, num)
  
  # Relax grouping: allow small gaps (max_gap)
  df <- df %>%
    arrange(letter, num) %>%
    group_by(letter) %>%
    mutate(group = cumsum(c(TRUE, diff(num) > max_gap))) %>%
    group_by(letter, group) %>%
    summarize(
      range = if (n() == 1) first(code) else paste0(first(code), "–", last(code)),
      sort_val = first(num) + utf8ToInt(first(letter)) * 1000,
      .groups = "drop"
    )
  
  list(
    range = paste(df$range, collapse = ", "),
    sort_key = min(df$sort_val)
  )
}

# Summarize and apply function
icd10_summary <- icd10categories %>%
  group_by(category_0321update) %>%
  summarize(
    `Number of ICD-10 codes` = n_distinct(icd10code),
    temp = list(condense_icd10_codes(icd10code, max_gap = 2)),
    .groups = "drop"
  ) %>%
  mutate(
    `ICD-10 code groups` = map_chr(temp, "range"),
    sort_key = map_dbl(temp, "sort_key")
  ) %>%
  rename(category = category_0321update) %>%
  arrange(sort_key)  # Alphabetical ICD-10 ordering

# Create GT table
gt_table <- icd10_summary %>%
  select(-sort_key, -temp) %>%
  gt() %>%
  tab_header(
    title = md("**Table 1. ICD-10 Code Groupings and Counts by Category**")
  ) %>%
  cols_width(
    `ICD-10 code groups` ~ px(500)
  ) %>%
  tab_options(table.width = pct(100))
gt_table
Table 1. ICD-10 Code Groupings and Counts by Category
category Number of ICD-10 codes ICD-10 code groups
infectious 199 A00–A99, B01–B99
cardiovascular 100 I00–I99
infectious respiratory 13 J00–J06, J09–J11, J20–J22
pneumonia 7 J12–J18
respiratory 78 J19, J23–J99
digestive 25 K20–K31, K50–K52, K55–K63, K92
injury 342 M48, M80, M84, M96–M97, R45, S00–S99, T00–T88, V00–V05, V09–V15, V18–V29, V40–V49, V53–V74, V79–V83, V86–V98, W00–W01, W04–W22, W25–W31, W34, W39–W40, W44–W46, W49–W57, W61, W67–W69, W74, W86–W88, X00–X02, X08, X14, X19, X30, X58, X74, X78–X80, X83, X95, X99, Y00, Y03–Y04, Y07–Y09, Y21, Y24, Y27–Y32, Y63, Y69–Y73, Y79, Y82–Y84, Y90–Y95, Y99
renal 29 N00–N28
other 75 O00–O04, O07–O16, O20–O36, O40–O48, O60–O77, O80–O82, O85–O94, O98–O99, Z31–Z33
covid 1 U07
# Save table
gtsave(gt_table, filename = "icd10_category_groupings.png")